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Qh The use of multigrid and related preconditioners with the finite element method is often 

limited by the difficulty of applying the algorithm effectively to a problem, especially when the 
domain has a complex shape or the mesh has adaptive refinement. We introduce a simplification 
of a general topologically-motivated mesh coarsening algorithm for use in creating hierarchies of 

i | meshes for geometric unstructured multigrid methods. The connections between the guarantees 

of this technique and the quality criteria necessary for multigrid methods for non-quasi-uniform 

h^- problems are noted. The implementation details, in particular those related to coarsening, 

remeshing, and interpolation, are discussed. Computational tests on pathological test cases 
from adaptive finite element methods show the performance of the technique. 

1 Introduction 

Multigrid methods enable the asymptotically optimal approximate numerical solution to a wide 
array of partial differential equations. Their major advantage is that they have, for a wide class of 
problems, 0(n) complexity for solving a linear system of n degrees of freedom (DoFs) to a specified 
accuracy. However, there are many stumbling blocks to their full implementation in conjunction 
with adaptive unstructured finite element methods. 
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Algebraic multigrid (AMG) methods and widely-available libraries implementing them have made 



great strides in allowing for the widespread use of multilevel methods in computational science 37 
This is due to their relatively black-box nature and ability to be used for a wide variety of common 
problems without much expert tuning. Interesting problems often must use meshes that have been 
adaptively refined in order to resolve numerical error or scale of interest without regard to multilevel 
discretization. In order to develop a fast multilevel method for such a problem, either an algebraic 
or geometric coarsening procedure must be employed. Practitioners faced with such problems are 
often immediately drawn to AMG because of the difficulties of geometric coarsening on nontrivial 
domains. 
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However, there are definite advantages to geometric multigrid. In fact, there have been several 
attempts to integrate some geometric intuition into algebraic methods 14p6p8 . For one, efficiency 



guarantees for AMG are based upon heuristic, computational observation, or theoretical bounds 
that are hard to satisfy computationally. Geometric multigrid methods, on the other hand, can be 
proven to have optimal convergence properties assuming certain conditions on the meshes, which 
will be elaborated upon in Section [L~2| 



An observation made by Adams [T] on the state of unstructured geometric multigrid methods was 



that while Zhang 39 40 provided guarantees of the optimality of geometric multigrid methods, 
they are not directly applicable to the methods of mesh coarsening and hierarchy construction as 
they are typically implemented. Our goal here is to use computational geometry methods to meet 
these quality criteria, to show a way of constructing a set of meshes in an optimal fashion that would 
satisfy these criteria, and to show that we have constructed a successful optimal-order multigrid 
method. 



1.1 Adaptive Finite Elements and Multigrid 

A fairly typical scheme for constructing a series of meshes for geometric multigrid methods involves 
uniform refinement from some coarsest starting mesh. However, the goals of refinement for the 
purpose of resolving discretization or modelling error and refinement for the purposes of creating 
a multigrid hierarchy can easily be at odds with each other if care is not taken. If a mesh has 
been refined to satisfy some analytical or adaptively determined error control, then it may have 
vastly different length scales existing in different parts of the mesh. It is because of this that the 
combination of refinement in order to resolve physics and refinement for the purpose of multigrid 
is fraught with peril. 

If one refines for the purpose of creating a multigrid hierarchy from an initial coarse mesh with 
some grading to handle the numerical properties of the problem, the refinement of the finest mesh 
will not reflect the error properly. Conversely if some stages of the refinement for the physics are 
used for the multigrid hierarchy, one may be unable to guarantee that the meshes would satisfy the 
quality and size constraints required by multigrid. Because of this it is very difficult to make the 
two concepts, adaptive refinement and geometric multigrid, go hand in hand as a single adaptive 
meshing strategy that satisfies the needs of both methods in the general case. 

In typical applications in computational science, one is often given a mesh that has already been 
adapted to the physics of some problem and asked to solve some some series of linear equations 
repeatedly on this static mesh. Two examples where this happens quite frequently are in optimiza- 
tion problems, and in the solution of nonlinear equations by methods requiring the equations to 
be linearized and solved repeatedly. In this case, the only available geometric multigrid methods 
are based upon coarsening, and huge advantages may be reaped from precomputation of a series of 
coarser spaces that effectively precondition the problem. 

The need for a mesh refined to resolve error can be demonstrated in some fairly straightforward 
cases. The need for a priori grading around a reentrant corner to resolve pollution effects [5] is 
a well-studied classical phenomenon [3 10 22 in adaptive finite element computation. Multigrid 



computations on reentrant corner problems have been analyzed on simple meshes in the two di- 



mensional case with shape-regular refinement 13 or structured grading 19 . A mesh arising from 
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these requirements has disproportionately many of its vertices concentrated around the reentrant 
corner, such as the refined mesh shown in Figure [I] Around a reentrant corner of mesh interior 
angle 9 r > it we can define a factor 

> 7r 
9 r 

Then, given some constants C a > and C < oo related to the maximum length scale in the mesh, 
a mesh that will optimally resolve the error induced by the reentrant corner, for h as the length 
scale of the cell containing any point a distance r away from the corner, will satisfy: 

Car 1 ' 1 " <h< Ctr 1 -" 

This problem will be considered further in Section [4j 





Figure 1: Quasi-uniform (right) and a priori (left) refined meshes: 500 vertices 



1.2 Unstructured Geometric Multigrid 

When the creation of a coarsened hierarchy is considered in the geometric sense, we mean that the 
problem is reduced to finding some series of meshes {M° ...M n } where M° is the initial mesh and 
and M" is the coarsest mesh. This hierarchy is constructed starting from M°. The differential 
operator is either rediscretized on each mesh in order to create the series of coarse problems, or 
interpolators between function spaces on the meshes are used to create approximate coarse versions 



of the original fine linear system, a process known as Galerkin multigrid 36 . The experiments in 
this paper are all done by rediscretization. 

One restriction that we will maintain throughout this paper is that the meshes are simplicial and 
node-nested. This means that the vertices of mesh M i+1 are some subset of the vertices of mesh 
M l . This is an assumption used by many mesh refinement and coarsening methods, as it is a 
reasonable and common assumption that the vertices of the initial fine mesh adequately capture 
the domain's geometry. This is in contrast to fully nested multigrid hierarchies, where each coarse 
cell is the union of several fine cells. A series of non-nested unstructured coarse meshes do not 
necessarily need to be node-nested either. For instance, the mesh hierarchy may be created by 
repeatedly remeshing some CAD or exact description of the domain at different resolutions. 

There is some justification for extending the classical multigrid convergence results to the case of 



non-nested meshes ISllQ 24 . Results of optimal multigrid convergence have been extended to non- 



quasi-uniform |40 and degenerate meshes 41 in two dimensions. Three dimensional non-nested 



multigrid has been proven to work in the quasi-uniform case 33 . 
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We will use the mesh quality conditions required for the non-quasi-uniform case 40 . While the 
example mesh hierarchies used in this paper might be quasi-uniform when considered alone, we may 
not assume quasi-uniformity independent of the size of the finest mesh in the asymptotic limit of re- 
finement. It should be noted that the non-quasi-uniform theory does not exist for three dimensions, 
but the quasi-uniform theory 



33 serves as a guide to how the method should perform. 



1.3 Mesh Constraints 

For any mesh M and any cell r in M let p T be r's incircle diameter and h T its longest edge, and 
define the aspect ratio as 

AR(r) = ^. (1) 

The quality of approximation 34 and matrix condition number [I] in FEM calculations depend 
strongly on AR(t). For some constant Car, we must require that the aspect ratio of any given cell 
in {M ...M n } satisfy 

AR(t) < C AR . (2) 

The other half of the criteria are the local comparability conditions. Assume that if we have two 
meshes, M k and M k ~ 1 , we define, for some cell r in M k 

S^{T:TeM M ,TnT^«} (3) 

which defines the set of cells in a mesh M fc_1 overlapping a single cell r in the next coarser mesh 
M k . We can state the local comparability conditions as 

sup {\S k \} < C for some C a > (4a) 

TGM fc 

sup { sup {t^-}} < Ci for some C; > 1 (4b) 



Eq. [4a] implies that each cell intersects a bounded number of cells in the next finer mesh, and Eq. 4b 



that the length scale differences of overlapping cells are bounded. We will also state here, for the 
sake of completeness, the assumption from the standard proofs of multigrid convergence necessary 
for the algorithmic efficiency of the method. Define 

\M\ = # of cells in mesh M 

Then, for some C m — 2 

\M k \>C m \M k+1 \. (5) 

The implication is that at each coarser mesh there must be a sufficient decrease in the number of 
cells. A geometric decrease in the number of cells over the course of coarsening is necessary in order 
to have an 0(n) method. 
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2 Mesh Coarsening 

For node-nested coarsening of an unstructured mesh, we separate the process of mesh coarsening 
into two stages: coarse vertex selection and remeshing. The problem of coarse vertex selection on a 
mesh M involves creating some subset V c of the vertices V of M that is, by some measure, coarse. 
Several techniques have been described for coarse vertex selection on unstructured simplicial meshes. 
The most widely used of these methods are Maximum Independent Set (MIS) algorithms [2| [TT[[T7] . 
These choose a coarse subset of the vertices of the mesh such that no two vertices in the subset are 
connected by a mesh edge. The resulting set of vertices is then remeshed. 

Given some mesh M, we define the graph Gm = (V, E), which consists of the edges E and vertices 
V of M. The MIS algorithms may then be described as finding some set of coarse vertices V^ IS 
such that 

{v a , Vb) E for all pairs of vertices (v a , Vb) € V C M/S . 

This may be implemented by a simple greedy algorithm where at each stage a vertex is selected 
for inclusion in V^. and its neighbors are removed from potential inclusion. There are a couple 
of issues with this method. The first being that there is no way to determine the size of the mesh 
resulting from coarsening. The other is that there are no real guarantees on the spatial distribution 
of vertices in the resulting point set. It has been shown that MIS methods may, for very simple 



example meshes, degrade the mesh quality quite rapidly 25 . Other methods for choosing the coarse 
vertex selection have been proposed to mitigate this shortcoming, often based upon some notion of 
local increase in length scale 130] . 



2.1 Function-Based Coarsening 

The coarse vertex selection method we will focus on for this was developed by Miller et al. [25] and 
is referred to as function-based coarsening. The method begins by defining some spacing function 
Sp(v) over the vertices v € V of M. Sp(v) is some approximation of the local feature size. In 
practice, the approximation of the local feature size based upon the nearest-mesh-neighbor distance 
at each vertex is a reasonable and efficiently computable spacing function. 

Define /3 as the multiplicative coarsening parameter. /3 can be considered the the minimum factor 
by which Sp(v) is increased at v by the coarsening process. Here we choose /3 > /3q where (3q = v2 
in 2D and /3o = \/3 in 3D by the simple fact that /3 = /3q + e for small e would reproduce the 
repeated structured coarsening of an isotropic, structured, shape-regular mesh where the length 
scale is increased by two in every direction. /3 may also be tuned, changing the size of the set 
of resulting coarse vertices in a problem-specific fashion to account for mesh and function-space 
properties, such as polynomial order. 

We say that a coarse subset of the vertices of M, V^ BC satisfies the spacing condition if 

/3(Sp(vi) + Sp(vj)) < dist(vi,Vj) for all pairs (vi,Vj) € V C FBC . (6) 

After determining V FBC , M c may be created by some remeshing process. A hierarchy of meshes 
may then be created by reapplying the coarsening procedure to each newly-minted coarse mesh 
with constant j3 and some recalculated Sp in turn to create a yet coarser mesh. This may be done 
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until the mesh fits some minimum size requirement, or until the creation of a desired number of 
meshes for the hierarchy. 



The original authors 25 had fast numerical methods in mind when proving quality bounds for 
function-based coarsening, and a number of the properties required of the meshes are spelled out 
in the original work. For one, the maximum aspect ratio of the resulting mesh hierarchy may be 
bounded by some constant, satisfying Eq. [2j We also have that, with hi(x) being the length-scale 
of the cell overlapping the point x for all x G Q: 

— hi + i(x) < hi(x) < Xhi+i(x) for some I > 0. 



Note that this condition implies the first of the local comparability conditions (Eq. 4a). By combin- 
ing the length scale bound with the aspect ratio bound, one may infer that only a certain number 
of fine cells may be in the neighborhood of a coarse cell, bounding the number of overlapping cells 



(Eq. 4b ) . Finally, for the coarsest mesh M n and some small constant b, 

\M n \<b. 

This implies that the coarsening procedure will be able to create a constant-sized coarse mesh. 
Because of these conditions and their parallels with the conditions on the multigrid hierarchy, this 
method for coarsening is particularly appealing. 

Parts of the function-based coarsening idea have been incorporated into other methods for mesh 



coarsening by Ollivier-Gooch 28 . Some similarities between the method we propose here and this 
previous work are that both use traversal of the mesh and local remeshing in order to coarsen 
meshes that exhibit complex features. This method constructs the conflict graph 

G C = {V,E C ) 

where 

Ec = {{vi,Vj) : P(Sp(Vi) + Sp{vj)) > dist(vi,Vj)}. 

This graph is then coarsened with an MIS approach as shown above. Note that in the limit of 
extreme grading or large j3 the Gq can grow to be of size 0(|y| 2 ). Our method avoids this by 
localizing and simplifying the notion of the spacing function and choice of vertex comparisons. We 
modify function based coarsening in a way that reliably guarantees optimal complexity irregardless 
of mesh non-quasi- uniformity as discussed in Section |1.1| without dependence on the mesh and the 
parameter ft. 



2.2 Graph Coarsening Algorithm 

We describe here a greedy algorithm for determining an approximation V ( F CA to the set of coarse 
vertices Vf BC satisfying a weakened notion of the spacing condition based upon Gm- Talmor 



38 proposed using shortest-weighted-path distance determined by traversing along mesh edges to 
accomplish this. Instead of doing using the shortest-weighted-path approach, we choose to use 
the edge connectivity to progressively transform Gm into some final G Mc , which approximates the 
connectivity of the coarse mesh, M c . It should be noted that, beyond the initial Gm formed at the 
start of the algorithm, Gm does not necessarily satisfy the condition that its edges represent the 
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(c) Subalgorithm on v (d) Subalgorithm on v\ 



Figure 2: This shows coarsening of a patch of Gm- Sp is defined, multiplied by /3, and the 
subalgorithm is invoked for two vertices, v and v\ shown in Figs. |2(c) and |2(d)| in dark grey with 
vertices violating the graph spacing condition shown in light gray. 
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edges of a valid simplicial mesh. Remeshing given M and V GCA as it is created in this section is 
discussed in Section [Ol 



We begin this process by modifying condition [6] to be 

P{Sp{v{) + Sp(v 2 )) < dist(v 1 ,v 2 ) for all {vi,v 2 ) G E of G M - (7) 

This restricts the spacing condition to take into account only distance between vertices connected 
by graph edges rather than all pairs of vertices. This restriction makes sense considering the 
calculation of the Sp(v) as nearest-neighbor distance is based upon adjacency of the mesh. This is 
an important simplification for complex domains, where the mesh may represent complex features 
of the domain such as holes, surfaces, and interfaces with different local characteristics on each side 
that may be hard to encode when the vertices are considered as a point cloud. In our model, the 
edges of the mesh are seen as encoding the topological connectivity of the domain. Spacing only 
applies to topologically connected areas of the mesh as encoded in Gm- 

Define 

F(y) = unknown, included, or excluded for alk? G V. 

The algorithm starts with F(y) = unknown for all v G V. Vertices that must be included for some 
reasons such as the boundary issues described in Section [2.4| are set to be included automatically 
and visited first. The remaining v £ V are then visited in some arbitrary order. If F(y) — excluded, 
the iteration skips the vertex. If F(v) = unknown, F(v) is changed to included. 

When F(v) is set to included, a subalgorithm is invoked that, given Gm and some particular vertex 
v, transforms Gm to some intermediate G M . This subalgorithm corresponds to coarsening the area 
around v until the spacing function is satisfied for all Nq v (v). Each w G Nq u (v) are tested to see 
if the edge (v, w) violates Eq. [7] In the case that the condition is violated and F(w) = unknown, 
w is removed from Gm and F(w) is changed to excluded. A removed w's neighbors in Gm are 
added to Nq m (v) by adding edges in Gm between v and all u G Ng u {w) if u £ Ng m (v) already. 
This may also be considered as an edge contraction from w to v. 

The outcome of this subalgorithm, G M , has that all w G Ng* iy) such that F(w) = unknown obey 
Eq. W\ There is the possibility that for some V\ G Ng^(v), F(vi) = included and Eq. [7] is not 
satisfied. This may arise due to the necessary inclusion of boundary vertices due to conditions 
described in Section [274} After the subalgorithm terminates, one is left with G V M . The algorithm 
then visits some v\ (Figure 



2(d) ), which, if included, will create some G M by modification of G M 



Once every vertex is visited, the whole algorithm terminates. At this point we may define 

v gca = | y . v e v,F(v) = included}. 

Despite the fact that we are considering coarsening and remeshing as a preprocessing step, we still 
have the goal that the algorithm should be 0(n) with n as the number of DoFs in the discretized 
system. For piecewise linear finite elements, the number of DoFs in the system is exactly the 
number of vertices, |V|, in M. For all reasonably well-behaved meshes, we can additionally assume 
that \E\ = 0(\V\) and \M\ = 0(\V\). This implies that n = 0(\V\). The complexity of a single 
invocation of the subalgorithm may be 0(|V| + \E\) if some large fraction of the mesh vertices are 
contracted to some v in one step. However, the aggregate number of operations taken to reach 
G Mc must be the order of the size of the original graph Gm ■ This aggregate complexity bound is 
independent of the order in which v are visited. While here we keep the assumed order in which 
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the vertices are visited arbitrary, we see in Section 2.4 that spccifing some restrictions on the order 
may be used to preserve mesh features. 

Theorem 2.1. Given a graph derived from a mesh Gm, the graph coarsening algorithm will create 
V C GCA in 0(\V\ + \E\) time. 

Proof. The fact that the complexity of the algorithm depends on \V\ at least linearly is apparent 
from the fact that every vertex is visited once. In order to show that the entire algorithm is 
0(|^| + \E\) we must establish that each edge e is visited at most a constant number of times 
independent of the size of the graph. 

As vertex v is visited and F(v) is set to included, the subalgorithm inspects each e = {v,n) for 
n e Ng m (v) to see if they satisfy the spacing condition with v. e is either deleted from Gm if n 
and v do not satisfy the spacing condition, or left in place if the spacing condition is satisfied or if 
F(n) = included. 

Edges that are deleted are not ever considered again. Therefore, we must focus on edges that survive 
comparison. Suppose that an edge e = (vo,i>i) is considered a second time by the subalgorithm at 
the visit to vertex vq. As this is the second visit to consider e, F(vi) must be included as in the 
first consideration of e necessarily happened during the visit to v\. As both endpoints of e have 
now been visited, there is no way e may be considered a third time. As each vertex is visited once 
and the distance across each edge in Gm is considered no more than twice, the algorithm runs in 
Q(\V\ + \E\) time. □ 



2.3 Remeshing 

In the node-nested case, the problem of remeshing may be treated as a procedure taking a mesh M 
and some coarse subset of M's vertices V c and returning a new mesh M c that approximates the 
domain spanned by M but only includes vertices V c . Traditional constrained remeshing methods 
using standard meshing technology require specialized information about the domain, including 
the location and topology of holes in the domain, that is not typically provided with a pre-graded 
fine mesh. This inhibits automation since a great deal of expert input is necessary to preserve the 
boundary during coarsening. Instead, we choose to locally remove the set of vertices V \V C from 
M one at a time to create M°. 

In the 2D case, we remove a vertices from M in-place by a simple series of geometric predicates 
and edge flips. This procedure is completely localized to the neighborhood around the vertex being 
removed. This is done in a way that retains the Delaunay condition 27 if the initial mesh was 
Delaunay. Similar but much more complex techniques exist in 3D [21 , but are difficult to make 
computationally robust enough for the requirements of multigrid methods, where the creation of a 
single bad cell may severely impact the strength of the resulting multigrid preconditioner. In this 
work, we relax the Delaunay condition and instead focus on mesh quality of the sort required by 



the multigrid method as stated in Section 1.2 



An extremely simple method of removing a vertex v from the mesh M is by quality-conserving edge 
contraction. Similar approaches have been proposed for this problem previously [28] . This may be 
stated as choosing n* e N(v) such that for the set of cells C n that would be created during the 
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contraction, 

n* = argrnin(rnax AR(c)) 



subject to the constraint 



n£iV(») ceC " 



max AR{c) < Car- 



Each resulting c must also be properly oriented in space. 

The computational work done to remove a single vertex is 0(|C„||iV(i>)|) as each potential configu- 
ration C n has that, when compared to the set of cells adjacent to v, C v , |C„| < \C V \ because at least 
two cells adjacent to both n and v must be removed by each contraction. For each n € N(v) one 
must check the aspect ratio and orientation of each cell in the associated C n . Assuming bounded 
mesh degree, the removal of 0(|y|) vertices would take 0(| V|) time. 

Note that we are directly enforcing the maximum aspect ratio that a given mesh modification may 
create, which makes this akin to best effort remeshing. If it is impossible to satisfy the quality 
conditions required of the mesh while removing a vertex, the vertex is left in place, at least for 
the time being. Such vertices are revisited for removal if their link is modified by the removal 
of one of their neighbors, a process that may only happen 0(N(v)) times, bounding the number 
of attempts that may be made to remove a given vertex. Due to this, we cannot guarantee that 
this 3D remeshing finds a tetrahedralization including just the coarse vertices. However, we note 
that our experience with this method is that the number of vertices that must be left in place is 
significantly less than the typical number of Steiner points added through direct Delaunay remeshing 



using Tetgen 35 , and that this method produces tetrahedralizations entirely suitable for use with 
multigrid methods. The results of a study of the quality of meshes created by this method is shown 
in Section |4~T1 



2.4 Mesh Boundaries 

Preservation of the general shape of the domain is important to the performance of the multigrid 
method as the coarse problem should be an approximation of the fine problem. In the worst 
case a coarser mesh could become topologically different from the finer mesh leading to very slow 
convergence or even divergence of the iterative solution. If this is to be an automatic procedure, 
then some basic criteria for the shapes of the sequence of meshes must be enforced. Therefore, the 
vertex selection and remeshing algorithms must be slightly modified in order to take into account 
the mesh boundaries. 

First, we must define the features, in 2D and 3D, which require preservation. We choose these 
features by use of appropriate-dimensional notions of the curvature of the boundary. Techniques 



like these are widely used in computational geometry and computer graphics 15 31 . In 2D, the 
features we choose to explicitly preserve are the corners of the mesh. We consider any boundary 
vertex with an angle differing from a straight line more than Cjc to be a corner. For the sake of 
this work we assume Ctc = § ■ 

In 3D, the discrete curvature at a boundary vertex is computed as the difference between 27r and 
the sum of the boundary facet angles tangent to that vertex. Vertices where the absolute value 
of the discrete curvature is greater than Ck. are forced to be included. High-dihedral-angle edges, 
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Figure 3: Interior and boundary coarsening stages. Dark gray vertices are included automatically; 
Light gray circles violate the spacing condition. 



corresponding to ridges on the surface of the domain, must also be preserved. We consider any 
boundary edge with the absolute value of the dihedral angle greater than Cjc to be a boundary 
ridge. 

Our approach to protecting the boundary during the coarse vertex selection procedure is to separate 
vertex selection into two stages, one for the interior, and one for the boundary. In the interior stage, 
the boundary vertices are marked as included automatically, and therefore any interior vertices 
violating the spacing condition with respect to a boundary vertex are removed from Gm (Figure 
3(a)[ ) . All boundary vertices now respect the spacing condition with respect to all remaining interior 



vertices, making the second stage, boundary vertex selection, entirely independent of the previous 
coarsening of the interior. The boundary vertex selection procedure then operates independently of 



this, producing the completely coarsened graph G A ~ f (Figure 3(c)). In 3D this process is repeated 
once again. During the boundary coarsening procedure, vertices lying on edges that have been 
identified as ridges are automatically included. The ridge vertices are then coarsened in turn. 
Corner vertices are automatically included during all stages of vertex selection, as shown in Figure 
|3(b)| For remeshing, the only modification necessary for boundary preservation is that the removal 
of a boundary or ridge vertex may only be done by contraction along an edge lying on the boundary 
or ridge. This simple modification works especially well for meshes that have a lot of flat planes on 
their surfaces with sharp ridges separating them. Without this procedure, ridges may appear torn 
or flipped in the coarser meshes, changing the shape of the domain significantly. As the hierarchy 
coarse meshes are created, new edges or vertices may become marked as corners or as belonging to 
ridges. 



3 Interpolation Operators 



We have shown that local traversal is a powerful tool for the construction of node-nested coarse 
meshes. Traversal of the node-nested meshes may also be used to efficiently construct the interpo- 
lation operators /; for I € [0, n— 1]. /; is defined as the interpolation operator between the finite 
element function spaces and Vf, which are defined over the coarse mesh M c and fine mesh 
respectively for the adjacent / = / and c = I + 1. We also assume that the domain £1 spanned by 
all M is connected. For the purposes of this paper, we assume that for each DoF i, associated with 
basis function M in V ft in the system has some nodal point X{ £ M. d and that a function / may be 
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projected into V£ by constructing the coefficients 

fi = f{xi) ■ rp{(xi). 

Given this assumption, we may construct the prolongation operator from functions in Vfi to V/ by 
associating each i with a cell t£ in M c and constructing, for ip{ and basis functions ip? supported 
on r fc c : 

One can, for the sake of simplicity, associate every ipf with a single cell r/ . The choice of cell is 
arbitrary between the cells the unknown is supported on. This makes the problem equivalent to 
finding, for some set of x € ft, T c (x), defined as the cell in M c that x is located in. We also define 
T c (x), which is some cell in M c that is, in some sense, nearby x. In addition to the points Xi 
associated with all ipf , we have the midpoint x s of every cell r/ G . 

We define a two- level nested-loop traversal algorithm. The outer loop locates T c (x T f) for each 
G , and the inner loop determines T c (xi) for all Xi associated with tp{ supported on t( . 
Once T c (xi) is resolved for all i, one may construct the nodal interpolation operator. 

The outer loop consists of breadth-first searches (BFSes) on the graph of the neighbor relations of 
the cells of M c in order to determine T c (x Tk ) for each £ . This is implemented by a simple 
FIFO queue in which the neighbors of a given cell are pushed on the back of the queue when it is 
visited. Enqueued and visited cells are marked as such and not enqueued again. We say that two 
cells Tq and rf are neighbors if they share a vertex, edge, or facet. The BFS to locate starts with 
the cell f c {x Tk ). 

As T c (x r f) is established for each cell t[, the inner loop is invoked. This loop consists of a BFS for 
each ip{ associated with t[ and determines T c (xi) by BFS starting from T c (xi) = T c (x r f). 

The last ingredient in the algorithm is how to determine T c (x /). One simple way of doing this is 
setting T c (x /) = T c (x s ) for any cell that is neighboring t[. This notion of locality may be 

k m 

exploited by setting the values of T c (t„) for all neighbors r„ of a cell r upon determining T c (t). 
When T c (t^) for any € is determined, the connectivity of M c and may be effectively 
traversed in tandem in order to extend the search over the whole meshed domain. 

This process may be both started and made more efficient by exploiting the node-nested properties 
of the meshes. We have that meshes and M c will share some set of vertices V^'^ due to our 
node-nested coarsening procedure. Define r/ to be a cell in mesh that is adjacent to some 
v G V"W' C ) and to be an arbitrary cell in M c adjacent to that same v. Then, one may initialize 
f c (x T f) to be t c v . 

An obervation about the complexity of this algorithm may be established assuming that the condi- 
tions in Section [l.2| are satisfied by the hierarchy of meshes. It should be obvious that the complexity 
is going to be bounded by the total number of BFSes done during the course of the algorithm mul- 
tiplied by the worst-case complexity of an invocation of the BFS. It's easy to see that for \M l \ cells 
in the fine mesh and n unknowns, the geometric search must be done \M \ + n times. 

In order to bound the complexity of a given BFS to a constant, we must show how many steps 
of traversal one must search outwards in order to locate some T C (W) given T c (t^). This may be 
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accomplished by showing that the length scale the search has to cover on the fine mesh may only 
contain a given number of coarse mesh cells. We know that the dist(x j , x T j ) is going to be less 
than some maximum search radius r T s = h T j + h / given that they arc adjacent. Wc also may put 

a lower limit on the minimum length scale of cells in M° that overlap and Tq by the constant in 
Eq.~ 



4b 



and by the aspect ratio limit in Eq. [2J as h^ in — , . 



_ c Q h Tf 



This gives us that the number of cells that may fit between T c (r^) and T c {t^) is at most r- 

which is independent of overall mesh size. Note that the distance between the center of a cell r 
and some nodal points Xi of i^i supported on r is always less than or equal to h T . This extends this 
constant-size traversal bound to the inner loop also, making the entire algorithm run in 0(|M^| +n) 
time. In practice on isotropic meshes T{t^) and T(t*) are almost always in the same cell or one 
cell over, meaning that the topological search is an efficient method for building the interpolation 
operators. 

An extension of the algorithm that has proven necessary on complex meshes is to allow for inter- 
polation to Xi not located within a cell of the mesh. For example, the Pacman mesh (Fig. [4J, when 
coarsened, will have vertices that lie outside the coarser mesh on the round section. In order to do 
this, the outer loop BFSes are replaced by a more complex arrangement where the BFSes search 
for the nearest cell rather than an exact location. This nearest cell then becomes T(x T ). The inner 
loop is replaced by a similar procedure, which then projects any xi that could not be exactly located 
to some Xi on the surface of the nearest cell and modifies the basis function projection to use <f) C j{xi) 
instead of (j)'j(xi). 

This procedure is easily extended to discretizations consisting of non-nodal DoFs by locating sets 
of quadrature points associated with the fine mesh cells instead of DoF nodes. This also allows for 
the construction of pseudointerpolants satisfying smoothness assumptions, such as the Scott-Zhang 
type interpolants (32). 



4 Model Problems 

We test the multigrid method using standard isotropic non-quasi-uniform examples from the liter- 
ature. The domains used are the Pacman in 2D and the Fichera corner in 3D. The problems and 
boundary conditions computed are designed to induce a corner singularity that makes refinement 
such as that described in Section [O] necessary 

The standard model problem for the Pacman domain is [3] designed to have a pollution effect at 
the corner singularity. The Pacman domain is a unit circle with a tenth removed. The reentrant 
angle is therefore 8 r = ~ radians. The associated differential equation is 

-Ait = in fl, 
u(r, 9) = ri sin on 

For the Fichera corner model problem, we separate the boundary of the domain into the reentrant 
surface dfli and the outer surface 9£! D . The Fichera corner mesh consists of a twice-unit cube 
centered at the origin with a single octant removed. The reentrant angle is 9 r = ^ radians 
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along the reentrant edges. A standard model problem with both Dirichlet and Neumann boundary 
conditions used for the Fichera corner domain is |29| 

in O, 

g on dfl , 

on <9fli. 

9xy + 9yz + Qzy i 

cos(#). 

where r and 9 consist of the radius and angle when restricted to the a;,-, Xj plane. Unfortunately, 
this problem does not have an analytical exact solution. Our goal is to show that we can effectively 
coarsen the a priori grading required to resolve the problem, so we will not be concerned with the 
solution accuracy as a measure of performance. 



-Au 
du 
dn 

U ■ 



Where 



9x iXj (r, 0) 



4.1 Mesh Quality Study 

In order to motivate the use of the method with the multigrid solver, one must look if our imple- 
mentation of the method actually lives up to the mesh quality metrics stated as requirements for 
guaranteed multigrid performance. This study was done using large initial Pacman and Fichera cor- 
ner meshes coarsened using the same algorithm used for the multigrid studies in Section |4.2| 

The measurements of mesh quality we have taken correspond to the bounds on the meshes in the 
multigrid requirements. These are the number of vertices and cells in the mesh, the worst aspect 
ratio in each mesh, the maximum number of cells in a mesh overlapping any given cell in the next 
coarser mesh, and the maximum change in local length-scale between a mesh and the next coarser 
mesh at any point in the domain. The meshes are created by repeated application of the coarsening 
algorithm with (3 = 1.5 in 2D and f3 = 1.8 in 3D. 



Mi 




Figure 4: Mesh hierarchy produced by function-based coarsening of an initial pacman mesh. 
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Tabic 1: Hierarchy quality metrics - 2D 



Pacman Mesh, 


6 = 1.5 


level 


cells 


vertices 


max( AR(t)) 


sup {\S k T \} 

reM k 


max 


( hk-i{x)} 
V h k {x) J 





1 

2 
3 
4 
5 
6 


70740 
22650 
7562 
2422 
811 
257 
86 


35660 
11474 
3858 
1254 
428 
143 
52 


3.39 
7.64 
7.59 
4.63 
5.52 
5.94 
7.76 


14 
15 
15 
15 
15 
12 


9.48 
6.23 
4.19 
5.32 
8.51 
4.94 




Figure 5: hierarchy of Fichera meshes 



Table 2: Hierarchy quality metrics - 3D 







Fichera Mesh, {3 = 1. 


3, C AR = 60 




level 


cells 


vertices 


max(AR(r)) 


sup {\S k T \} 


/ h k -i(x)\ 
max - 








T<£M k 


xen \ h k {x) J 





373554 


72835 


4.41 






1 


49374 


9120 


59.3 


197 


6.63 


2 


11894 


2269 


59.2 


131 


6.70 


3 


3469 


693 


59.3 


94 


6.68 


4 


914 


208 


58.8 


121 


6.61 


5 


182 


52 


49.5 


94 


6.87 



In 2D (Table [T]) the aspect ratio of the resulting cells stays within acceptable limits during the 
entire process. The slight increase in the aspect ratio is expected for the coarsening of highly 
graded meshes, as the coarser versions must necessarily remove vertices lying between refined and 
coarse regions of the non-quasi-uniform mesh. However, the associated increase in aspect ratio is 
kept to a minimum by the enforcement of the Delaunay condition. 

In 3D (Table [2]), we see consistent decrease of the mesh size and increase in the length scale. The 
maximum aspect ratio stays around Car for most of the levels. Further work on our incredibly 
simple remeshing scheme should be able to improve this. However, we do not see successive degra- 
dation of the quality of the series of meshes as the coarsening progresses. We can assume that the 
quality constraints are reasonably met in both 2D and 3D. 
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Table 3: Multigrid Performance on the 2D (Pacman) and 3D (Fichera) problems. Levels is the 
number of mesh levels. MG is the number of multigrid V-cycles. ILU is ILU-preconditioned GMRES 
iterations. L2 Error is computed based upon difference from the exact solution for the Pacman test 
problem. 

(a) Pacman Problem Performance (b) Fichera Problem Performance 



DoFs 


levels 


MG 


ILU 


Li Error 




DoFs 


levels 


MG 


ILU 


762 


3 


6 


76 


9.897e-6 




2018 


2 


7 


46 


1534 


4 


6 


129 


3.949e-6 




3670 


2 


8 


60 


3874 


5 


6 


215 


9.807e-7 




9817 


3 


8 


112 


9348 


5 


7 


427 


2.796e-7 




27813 


4 


9 


207 


23851 


6 


7 


750 


1.032e-7 




58107 


5 


9 


385 


35660 


7 


7 


1127 


8.858e-8 




107858 


6 


9 


543 


78906 


7 


7 


1230 


3.802e-8 




153516 


6 


9 


440 


139157 


8 


9 


3262 


1.535e-8 




206497 


7 


9 


616 


175602 


8 


9 


3508 


3.440e-9 




274867 


8 


9 


652 



4.2 Multigrid Performance 

The performance of multigrid based upon the mesh series created by this algorithm using the two 
test examples was carried out using the DOLFIN [23] finite element software modified to use the 
PCMG multigrid preconditioner framework from PETSc [6j. The operators were discretized using 
piecewise linear triangular and tetrahedral finite elements. The resulting linear systems were then 
solved to a relative tolerance of 10~ 12 . The solvers used were ILU- preconditioned GMRES for the 
standard iterative case, shown in the ILU columns as a control. In the multigrid case we chose 
to use GMRES preconditioned with V-cycle multigrid. Three pre and post-smooths using ILU as 
the smoother were performed on all but the coarsest level, for which a direct LU solve was used. 
For this to make sense, the coarsest mesh must have a small, nearly constant number of vertices; 
a condition easily satisfied by the mesh creation routine. We coarsen until the coarsest mesh is 
around 200 vertices in 2D or 300 vertices in 3D. 

These experiments show that the convergence of the standard iterative methods becomes more and 
more arduous as the meshes become more and more singularly refined. The singularity in 2D is 
much greater, due to the sharper reentrant angle, than it is in 3D, so the more severe difference in 
performance between multigrid and the control method is to be expected. We see that the number 
of multigrid cycles levels out to a constant for both 2D (Table 3(a) ) and 3D (Table [3(b) \ . We also 



see that a steadily increasing number of multigrid levels are automatically generated and that the 
method continues to behave as expected as the coarsening procedure is continually applied. 



5 Discussion 

We have proposed and implemented a method for mesh coarsening that is well suited to the construc- 
tion of geometric unstructured multigrid on 2D and 3D meshes. The technique used to construct 
the mesh hierarchy is novel in that it combines the advantages of graph-based mesh coarsening 
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algorithms with those of techniques considering the local feature size of the mesh. This is done in 
a straightforward way that only considers the local connectivity of the mesh. 

The further extension of the method to a parallel environment is fairly straightforward. Incor- 
porating a separate parallel coarsening stage for vertices on the processor boundaries would be a 
natural extension of the method. Techniques with a localized coarser solve, such as a variant of the 
Bank-Hoist method [7], could be easily replicated in this framework as well. 

Meshing issues still exist, as the local remeshing can be quite costly in 3D. There are several ways 
the remeshing could be improved, both in terms of efficiency and mesh quality guarantees. It is an 
experimental observation that remeshing using meshing software with constrained interior vertices 
is often faster than vertex removal. A mixed approach, where the domain could be topologically 
traversed and divided into regions that could be remeshed separately, may be an interesting compro- 
mise. Both approaches to remeshing are also easily extendable to the creation of an all-encompassing 
parallel algorithm. 

Despite this cost, the big win is twofold: treating the construction of coarse spaces as a precom- 
putation step in some nonlinear or optimization solve, and using higher-order elements, where the 
topological description of the mesh is smaller compared to the number of DoFs. It's been shown 
to be quite important to preserve the order of the approximation space when using multigrid for 



higher-order problems 20 , and the techniques for efficiently building both pointwise and other 
types of interpolants shown here already work in a similar fashion for higher order FEM discretiza- 
tions. As we can easily control the rate at which the mesh is decimated, we can coarsen more 
aggressively for the high order methods, retaining the optimal complexity of the method. 

One final future direction is to look at how this techique applies to problems on anisotropic meshes 
and the associated theory. There are successful node-nested technques proposed for semicoarsening 



anisotropic meshes 26 , and the coarsening technique described here could be applied in interesting 
ways. The other related problem of interest involves multigrid methods for problems with highly 
variable coefficients. Often these methods are developed and implemented in some nested framework 



12 , but with our ability to build meshes and interpolation operators that preserve interfaces (due 
to their construction by traversal) we should be able to look at these problems in the generalized 
unstructured setting. 



References 

[1] Mark Adams, Evaluation of three unstructured multigrid methods on 3D finite element prob- 
lems in solid mechanics, Int. J. Numer. Meth. Engng, 55 (2000), p. 519534. 

[2] Mark Adams and James W. Demmel, Parallel multigrid solver for 3D unstructured finite 
element problems, in Supercomputing '99: Proceedings of the 1999 ACM/IEEE conference on 
Supercomputing (CDROM), New York, NY, USA, 1999, ACM, pp. 27+. 

[3] Thomas Apel and Frank Milde, Comparison of several mesh refinement strategies near 
edges, Communications In Numerical Methods In Engineering, 12 (1996), pp. 373-381. 

[4] I. Babuska and A. K. Aziz, On the angle condition in the finite element method, SIAM 
Journal on Numerical Analysis, 13 (1976), pp. 214-226. 



GRADED UNSTRUCTURED MG 



18 



[5] Ivo Babuska, R. B. Kellogg, and J. Pitkaranta, Direct and inverse error estimates for 
finite elements with mesh refinements, Numerische Mathematik, 33 (1979), pp. 447-471-471. 

[6] Satish Balay, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. 
Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang, PETSc Web 
page, 2010. http://www.mcs.anl.gov/petsc. 

[7] Randolph Bank, Some variants of the Bank-Hoist parallel adaptive meshing paradigm, Com- 
puting and Visualization in Science, 9 (2006), pp. 133-144. 

[8] M. L. Bittencourt, C. C. Douglas, and R. A. Feijoo, Non-nested and non- structured 
multigrid methods applied to elastic problems, Part II: The three-dimensional case. Submitted 
for publication, (1998). 

[9] , Adaptive non-nested multigrid methods, Engineering Computations, 19 (2002), pp. 158- 

176. 

[10] H. Blum AND R. RANNACHER, Extrapolation techniques for reducing the pollution effect of 
reentrant corners in the finite element method, Numerische Mathematik, 52 (1988), pp. 539- 
564. 

[11] T. Chan, B. Smith, and J. Zou, Multigrid and domain decomposition methods for unstruc- 
tured meshes, 1994. 

[12] L. Chen, M. Holst, J. Xu, and Y. Zhu, Local Multilevel Preconditioners for Elliptic 
Equations with Jump Coefficients on Bisection Grids, ArXiv e-prints, (2010). 

[13] L. Chen and C. Zhang, A coarsening algorithm on adaptive grids by newest vertex bisection 
and its applications, Journal of Computational Mathematics, 28 (2010), pp. 767-789. 

[14] Edmond Chow, An unstructured multigrid method based on geometric smoothness, American 
Journal of Mathematics, 65 (2001), pp. 197-215. 

[15] N. Dyn, K. Hormann, S.-J. Kim, and D. Levin, Optimizing 3D Triangulations Using 
Discrete Curvature Analysis, Vanderbilt University, Nashville, TN, USA, 2001, pp. 135-146. 

[16] D. Feuchter, I. Heppner, S. Sauter, and G. Wittum, Bridging the gap between ge- 
ometric and algebraic multigrid methods, Computing and visualization in science, 6 (2003), 
pp. 1-13. 

[17] Herve Guillard, Node-nested multi-grid method with Delaunay coarsening, Research Report 
RR-1898, INRIA, 1993. 

[18] Jim E. Jones and Panayot S. Vassilevski, AMGE Based on Element Agglomeration, 
SIAM Journal on Scientific Computing, 23 (2001), pp. 109 133. 

[19] M. Jung, Some multilevel methods on graded meshes, Journal of Computational and Applied 
Mathematics, 138 (2002), pp. 151-171. 

[20] Michael Koster and Stefan Turek, The influence of higher order FEM discretisations on 
multigrid convergence, Computational Methods in Applied Mechanics, 6 (2005), pp. 221-232. 

[21] Hugo Ledoux, Christopher M. Gold, and George Baciu, Flipping to robustly delete a 
vertex in a delaunay tetrahedralization., in ICCSA (1), 2005, pp. 737-747. 



GRADED UNSTRUCTURED MG 



19 



[22] XlAOHAl LlAO AND RiCARDO H. Nochetto, Local a posteriori error estimates and adaptive 
control of pollution effects, Numerical Methods for Partial Differential Equations, 19 (2003), 
pp. 421-442. 

[23] A. LOGG and G. N. Wells, DOLFIN: Automated finite element computing, ACM Transac- 
tions on Mathematical Software, 37 (2010), pp. 20:1-20:28. 

[24] R. LOHNER AND K. MORGAN, An unstructured multigrid method for elliptic problems, Inter- 
national Journal for Numerical Methods in Engineering, 24 (1987), pp. 101-115. 

[25] Gary L. Miller, Dafna Talmor, and Shang-Hua Teng, Optimal coarsening of unstruc- 
tured meshes, J. Algorithms, 31 (1999), pp. 29-65. 

[26] E. Morano, D. J. Mavriplis, and V. Venkatakrishnan, Coarsening strategies for un- 
structured multigrid techniques with application to anisotropic problems, SIAM J. Sci. Comput., 
20 (1998), pp. 393-415. 

[27] M. MOSTAFAVI, Delete and insert operations in Voronoi/Delaunay methods and applications, 
Computers & Geosciences, 29 (2003), pp. 523-530. 

[28] Carl Ollivier-Gooch, Coarsening unstructured meshes by edge contraction, Int. J. Numer. 
Mcth. Engng., 57 (2003), pp. 391-414. 

[29] W. Rachowicz, D. Pardo, and L. Demkowicz, Fully automatic hp-adaptivity in three 
dimensions, Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 4816- 
4842. 

[30] B. Rey, K. Mocellin, and L. Fourment, A node-nested Galerkin multigrid method for 
metal forging simulation, Computing and Visualization in Science, 11 (2008), pp. 17-25. 

[31] Szymon Rusinkiewicz, Estimating curvatures and their derivatives on triangle meshes, 
3D Data Processing Visualization and Transmission, International Symposium on, (2004), 
pp. 486-493. 

[32] L. Ridgway Scott and Shangyou Zhang, Finite element interpolation of nonsmooth func- 
tions satisfying boundary conditions, Mathematics of Computation, 54 (1990), pp. 483-493. 

[33] , Higher- dimensional nonnested multigrid methods, Mathematics of Computation, 58 

(1992), pp. 457- 466. 

[34] Jonathan Richard Shewchuk, What is a good linear element? - interpolation, conditioning, 
and quality measures, in In 11th International Meshing Roundtable, 2002, pp. 115-126. 

[35] Hang Si, TetGen website, 2007. http://tetgen.berlios.de/index.html. 

[36] Barry Smith, Petter Bjorstad, and William Gropp, Domain Decomposition: Parallel 
Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 
March 2004. 

[37] K. Stuben, A review of algebraic multigrid, J. Comput. Appl. Math., 128 (2001), pp. 281-309. 

[38] Dafna Talmor, Well-Spaced Points for Numerical Methods, PhD thesis, Carnegie Mellon 
University, Pittsburgh, Pennsylvania, 1997. 



GRADED UNSTRUCTURED MG 



20 



[39] Shangyou Zhang, Optimal-order nonnested multigrid methods for solving finite element equa- 
tions I: On quasi-uniform meshes, Mathematics of Computation, 55 (1990), pp. 23-36. 

[40] , Optimal-order nonnested multigrid methods for solving finite element equations II: On 

non- quasi-uniform meshes, Mathematics of Computation, 55 (1990), pp. 439-450. 

[41] , Optimal-order nonnested multigrid methods for solving finite element equations III: On 

degenerate meshes, Mathematics of computation, 64 (1995), pp. 23-49. 



